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Abstract 

We analyze the dynamics of a model of a nanobeam under compression. The model is a two mode 
truncation of the Euler-Bernoulli beam equation subject to compressive stress applied at both ends. 
We consider parameter regimes where the first mode is unstable and the second mode can be either 
stable or unstable, and the remaining modes (neglected) are always stable. Material parameters 
used correspond to a silicon nanobeam. The two mode model Hamiltonian is the sum of a (diagonal) 
kinetic energy term and a potential energy term. The form of the potential energy function suggests 
an analogy with isomerisation reactions in chemistry, where 'isomerisation' here corresponds to a 
transition between two stable beam configurations. We therefore study the dynamics of the buckled 
beam using the conceptual framework established for the theory of isomerisation reactions. When 
the second mode is stable the potential energy surface has an index one saddle and when the 
second mode is unstable the potential energy surface has an index two saddle and two index one 
saddles. Symmetry of the system allows us to readily construct a phase space dividing surface 
between the two "isomers" (buckled states); we rigorously prove that, in a specific energy range, 
it is a normally hyperbolic invariant manifold. The energy range is sufficiently wide that we can 
treat the effects of the index one and index two saddles on the isomerisation dynamics in a unified 
fashion. We have computed reactive fluxes, mean gap times and reactant phase space volumes 
for three stress values at several different energies. In all cases the phase space volume swept out 
by isomerizing trajectories is considerably less than the reactant density of states, proving that 
the dynamics is highly nonergodic. The associated gap time distributions consist of one or more 
'pulses' of trajectories. Computation of the reactive flux correlation function shows no sign of a 
plateau region; rather, the flux exhibits oscillatory decay, indicating that, for the 2-mode model in 
the physical regime considered, a rate constant for isomerization does not exist. 

PACS numbers: 05.45.-a, 62.25.Fg, 62.25.-g, 82.20.-w, 82.20.Pm, 82.20.Sb, 82.20.Db, 83.10.Ff, 89.90.+n 
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I. INTRODUCTION 



There is currently much interest in the mechanical properties of nanoscale objects such as 
rods and cantilevers^-. For example, changes in resonant frequencies with mass loading can 
enable sensitive detection of molecular species with given mass. The possible manifestation 
of quantum effects are also of great interest (see ref. for a recent review). The effects of 
mechanical stress on nanostructures, such as the buckling of nanobeams, has been studied 
both experimentally and theoretically^-. This latter work is directly related to the work in 
this paper. 

The system studied in the present paper is a nanobeam subject to compressive stress. 
According to standard continuum mechanics, the beam will buckle as the magnitude of the 
associated compressive strain increases. More specifically, following previous work- - — , we 
consider a regime in which the dynamics of the buckled beam is usefully described by a 
2-mode model obtained by truncation of the full dynamics. In the regime considered, one of 
the modes is always unstable, while the second mode can be either stable or unstable, and 
all remaining modes (which are not considered explicitly) are stable. 

We derive a 2 degree-of-freedom (DoF) Hamiltonian describing the 2-mode dynamics. The 
Hamiltonian describes a bistable (reactive) mode coupled to a transverse degree of freedom; 
the dynamical system thereby obtained has precisely the form of simple model potentials 
that have been used to describe isomerization reactions in chemistry^— . However, for the 
nanobeam problem, the simple form of the potential is a rigorous consequence of the 2- 
mode truncation of the dynamics rather than an approximation to an unknown molecular 
potential energy function. 

We can therefore investigate the dynamics of the buckled nanobeam using the conceptual 
framework established for the theory of isomerization reactions. A number of approaches 
to the study of isomerization dynamics have been developed in the chemical literature (see, 



for example, refs ll6l- ll8U20ll2ll ). Some basic relevant questions are the following: can 'iso- 



merization' of the nanobeam be characterized by a ratel If so, can the rate coefficient 
be predicted using standard (so-called statistical) theories^^i, such as transition state 
theory^^^^^? If not, what are the dynamical properties of the system that lead to 
non-statistical reaction dynamics? 

Chakraborty et al. have previously applied harmonic transition state theory (TST) to 
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predict isomerization rates for various nanobeams under compressive stress^ - — . The rates 
obtained using basic harmonic TST are proportional to the inverse of the curvature of the 
potential transverse to the reaction coordinate, and so diverge at the strain value for which 
the second mode passes from being stable to unstable (and the associated saddle point 
passes from being index 1 to index ). Chakraborty has also applied a quantum version 
of harmonic TST to the nanobeam isomerization kinetics^— . 

In the present paper we examine the dynamics of a compressed buckled Silicon (Si) 
nanobeam from the perspective of reaction rate theory. The potential function which 
emerges from a standard modal analysis of the transverse beam displacements has a high 
degree of symmetry, and the appropriate dividing surface for computation of reactive flux is 
then in fact completely determined by symmetry. Moreover, we are able to rigorously prove 
that the dividing surface is a normally hyperbolic invariant manifold in a specific energy 
range. The energy range is sufficient for treating the effects of the index one and index two 
saddles on the isomerisation dynamics in a unified fashion. 

We compute both the distribution of gap times and the reactive flux at a number of 
energies in a physically relevant range. Our results show that, in this regime for the Si 
nanobeam considered: (i) the isomerization dynamics is extremely regular and nonergodic, 
and (ii) a rate constant for isomerization does not exist. Rather than exhibiting a rapid 
drop to a 'plateau' value followed by slow exponential decay^ 1 ^^, the reactive flux shows 
damped oscillatory decay. The system considered is in a regime (T > 100 K) where quantum 
effects are likely to be negligible. 

The structure of the paper is as follows. In Sec. [Ill we derive the equations of motion for the 
nanobeam. In Sec. Ill II we discuss in further detail the two degree-of-freedom Hamiltonian 
system obtained for the 2-mode truncation of the beam dynamics. Particular attention 
is given to the phase space geometry, specifically, the existence of a normally hyperbolic 
invariant manifold (NHIM)^ in the system phase space. In Sec. [IV] we review some concepts 
from reaction rate theory: phase space dividing surfaces and volumes, gap times and reactive 
fluxes. In Sec. [V] we discuss the physical parameter values and energy scales appropriate 
to our calculations. Sec. IVII presents the results of our numerical calculations, and Sec. IVHI 
concludes. In Appendix [A] we apply the concept of exponential dichotomies^ to provide a 
proof of the existence of a NHIM in the phase space of the truncated nanobeam problem. 
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II. EQUATIONS OF MOTION FOR THE NANOBEAM 



In this section we derive a 2-mode (classical) model for an Euler-Bernoulli beam subject 
to compressive stress applied at both ends. The derivation of the Euler-Bernoulli equations 
can be found in many textbooks on continuum mechanics (e.g. ref. |4|). Our derivation 



and notation closely follows that of refs l9Hlll (see also refs . A useful discussion of the 
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concept of stress in a quantum mechanical system is given in ref 

We consider an Euler-Bernoulli model of a beam of length L having width w and thick- 
ness d, where L ^> w > d. The requirement w > d allows us to assume that transverse 
displacements, y(x,t), occur only in the d direction. The linear modulus F (dimensions of 
energy per unit length) is related to the elastic modulus Q by F = Qwd. For a beam of 
rectangular cross-section the bending moment is given by k = C and /i = ^ denotes the 
mass per unit length. The length of the uncompressed rod is denoted by L . 

Constant compressive stress is applied to both ends of the beam, reducing the horizontal 
distance between the two endpoints of the beam to L < L . The strain e is e 



■^a, and is 



negative for compression. The compression causes a contribution to the potential energy of 
the beam due to bending in the d direction (the first term in (12. ip ) and elasticity (the last 
three terms in (12. ip ). where the potential energy has the form: 



V[y(x,t)} 



The kinetic energy is: 



dx {F K 2 (y") 2 + Fe(yf) 



8Lr 



L n 2 

l\2 



dx{y 



FL 



V 



T[y{x,t)\ = ^J dx ^V 
Forming the Lagrangian in the usual manner: 



L[y(x, t), y{x, t)] = T[y(x, t)} - V[y{x, £)] 



Lagrange's equations of motion are given by: 

d 5L 5L 
dt Sy 5y 



Using d2HD and (12^2]) . we have: 

ST[y(x,t)] 
Sy(x,t)] 



SV\y(x,t)] 
5y(x,t)\ 



Fk 2 V W 



Fey" 



2L 



dx(y'f y" 



(2.1) 
(2.2) 

(2.3) 

(2.4) 

(2.5a) 
(2.5b) 
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Using these expressions, together with (12. 3p and ( 12. 4p . gives the equations of motion: 



Fey" 



2Lr 



o \Jo 



dx{y') 2 y" 



(2.6) 



and the boundary conditions are chosen to be: 



y(0,t) = y(L,t) = 0, 
y"(0,t) = y"(L,t) = 0, 



(2.7a) 
(2.7b) 



which are referred to in the literature as hinged boundary conditions. 

One easily sees by inspecting (12. 6p that y(x,t) = is a solution. Linearizing ( 12. 6 p about 
this solution gives: 

jjjj = -Fn 2 y [i] + Fey". (2.8) 

We seek the normal modes (eigenfunctions) of these linearized equations by assuming a 
solution of the form: 

y n (x,t)=y n (x)e itJnt , n = 1,2,3,... (2.9) 

where 



y n (x) 



2 

— sin 



nirx 



L L L . 

Note that the y n (x) satisfy the normalization condition: 

y n (x)y m (x)dx = 5 n 
Substituting ([M} into §£E§ gives: 



(2.10) 



(2.11) 



IMJ 2 n y n (x) = F K 2 yW(x) - Fey'^ 



x) 



(2.12) 



We substitute (12.101) into (I2.12p . and after some algebra we obtain: 



with 



2 K F 
K 2r K 2 



(2.13) 



(2.14a) 
(2.14b) 
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Note that the quantity e is approximately equal to the critical value of the strain e c , obtained 
by solving the implicit equation (12.14b|) . but has a weak dependence on e. 

From the form of (12.91) and (I2.13p . we see that the mode y n (x,t) is linearly stable (resp., 
unstable) provided n 2 — I > (resp., n 2 — | < ). 

We will examine the situation where: 



1 - - < 0, 

e 

4 — ~ > or 4 — - < 0, 

e e 

n 2 - - > 0, n > 3. 



(2.15a) 
(2.15b) 
(2.15c) 



In other words, we will consider the cases where the first mode is always unstable, the second 
mode can be either stable or unstable, and modes n > 3 are all stable. 
Assuming that the solution of (I2.6P has the form: 

T 00 



y(x,t) 



sin 



nirx 



(2.16) 



n=l 



we substitute into ( 12.61) to obtain an infinite set of ordinary differential equations for the 
time evolution of the modal amplitudes, A n (t). However, we will simplify the problem by 
only considering the evolution of the first two modes: 



y(x,t) 



Ai(£)sin 



TTX 

lT 



A 2 (t) sin 



2nx 



(2.17) 



In this case one obtains a two-degree-of-freedom system for the evolution of the modal 
amplitudes Ai(t) and A 2 (t). Defining momentum variables Pi = fiAi, i = 1,2, the time 
evolution of the amplitudes A{ is described by a two degree-of-freedom Hamiltonian system, 
with Hamiltonian: 

H 



where the potential energy has the form: 



V(A 1 ,A 2 



F7T 2 (e 



-A\ 



2Fn 2 (e - 4e) 



AU 



Fir 4 



{A( + AA 2 2 



2\2 



(2.18) 



(2.19) 



2L 2 1 L 2 1 8L 4 L 

It is natural to ask how well the two-mode truncation described by the two degree- 
of-freedom Hamiltonian system defined by Hamiltonian ( I2.18P models the solution of the 
full partial differential equation describing the Euler-Bernoulli beam given in ( 12. 6p . It is a 
standard engineering approximation to approximate the full solution of a partial differential 



equation by considering a truncated modal expansion of eigenfunctions obtained from the 
linearized equations about an equilibrium state. The reasoning is that the evolution near 
the equilibrium solution is dominated by the evolution of the unstable modes. In some 
cases this can be rigorously proven using center manifold or inertial manifold techniques 
(see, e.g., refs l3llJ32l ). A seminal example of this approach that played a fundamental role 
in the development of applied dynamical systems theory was the work of Holmes, Marsden, 
and Moon on the dynamics of a buckled beam subject to periodic (temporal) forcing^—. 
Initially, a combination of experimental and theoretical work showed that the experimentally 
observed chaotic behavior was captured by the evolution of the one unstable mode, subject 
to forcing22i2i. Later, it was rigorously shown^ that this single mode truncation captured the 
dynamics of the full partial differential equation governing the beam (near the instability). 
In this paper we will not be concerned with these issues. Rather, we take as the starting 
point of our analysis the two degree-of-freedom Hamiltonian system governing the two mode 
truncation of the Euler-beam equation given in f ]2.6j) . 
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III. TWO-MODE TRUNCATION: HAMILTONIAN AND PHASE SPACE GEOM- 
ETRY 



We begin by non-dimensionalizing the two degree-of-freedom Hamiltonian system of eq. 
f)2.18p . Defining the dimensionless amplitudes: 



7T 



A 2 = ^A 2 



and substituting these expressions into the potential function (I2.19P gives: 



6 - e)A\ + (e - 4e)4^ + ±(Al + AA\f 



V(A 1 ,A 2 ) = FL 

= FL V(A 1 ,A 2 ). 
Defining associated momenta {pk} conjugate to the {A^} via 



Pi = Pi, 



71 



Pi = P2 



71 



and substituting into (12.181) gives the scaled Hamiltonian: 

H 

where 



JL = eL + & + v(a 1 ,a 2 ) 



2FL 2 L 2 , 



7T- 



Rescaling the momenta 



Pi ~^ nl/2 : 



% = 1,2 



we obtain the following Hamiltonian: 



with 



and 



H = § + § + V(Ai,A 2 ) 



2\2 



a = e — e 
/3 = e-4e. 



(3.1a 
(3.1b 

(3.2a 
(3.2b 

(3.3a 
(3.3b 

(3.4 

(3.5 
(3.6 

(3.7 
(3.8 

(3.9a 
(3.9b 
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The corresponding Hamiltonian equations of motion (for suitably rescaled time) are then: 

t dH , 

Ai = — =px, 3.10a 
dpi 

dH 

Pi = = -2Aj (a + A\ + 4A 2 2 ) , (3.10b) 
dH 

A 2 = — =p 2 , (3.10c) 
dH 

p 2 = -TTT- = ~8A 2 (P + A$ + AA 2 2 ) . (3.10d) 



A. Equilibria and their stability 

From the (algebraically) simple form of Hamilton's equations given in eq. (I3.10p it is 
straightforward to compute the equilibria, determine their linearized stability properties 
(i.e. compute the eigenvalues of the matrix associated with the linearization of Hamilton's 
equations about the equilibrium point), and compute the (total) energy of the equilibrium 
point. These properties are summarized in Table [H 

For a < 0, (3 > there are only three equilibrium points. The origin is an index one 
saddle point and the remaining two equilibria have two pairs of purely imaginary eigenvalues 
(minima of the potential ( I3.8P ). For a < 0, (3 < 0, |a| > \f3\, the origin is an index 
two saddle, phase space points (Ai,pi, A 2 ,p2) = (±\/— a, 0, 0, 0) have two pairs of purely 
imaginary eigenvalues and correspond to minima of the potential (I3.8p . and (^4i,pi, A 2 ,p2) = 
^0, 0, ±^^-, Oj are index one saddles. 

B. Invariant planes and the existence of a Normally Hyperbolic Invariant Manifold 

It can be seen by inspection of f)3. IQj) that, if we set A 2 = p 2 = (resp., Ax — pi — 0), 
then A 2 = p 2 = (resp., Ai = pi — 0). It then follows that the two planes: 

n x = {{Ai,Pi,A 2 ,p 2 )\A 2 = p 2 = 0} , (3.11a) 

n 2 = {{A 1 ,p 1 ,A 2 ,p 2 )\A 1 =p 1 = o} (3.iib) 

are each invariant with respect to the dynamics generated by f)3.10p . The dynamics on Hi 
is given by the Hamiltonian system defined by the Hamiltonian 

Hi^ + aA\+ l -A{ (3.12a) 
10 



and the dynamics on II2 is given by the Hamiltonian system defined by the Hamiltonian 



Hence, the dynamics on each plane is integrable. However, the dynamics on each plane is 
not isoenergetic. The three dimensional energy surface intersects a two dimensional plane in 
the four dimensional phase space in a one dimensional set, i.e. a trajectory of the one degree- 
of-freedom Hamiltonian system defined by H\ (for intersections with Hi) or a trajectory of 
the one degree-of- freedom Hamiltonian system defined by H 2 (for intersections with n 2 ). 

We now want to determine conditions under which some portion of n 2 is a normally 
hyperbolic invariant manifold (NHIM). Roughly speaking, NHIMs have saddle-like stability 
properties in directions transverse to the invariant manifolds^. In recent years NHIMs have 
been shown to be a significant phase space structure related to reaction dynamics. For 
example, they play the key role in the construction of a phase space dividing surface having 
the no-recrossing property^ 1 ^ and minimal flux—. They have also been shown to be central 
to Thiele's theory^ of reaction dynamics in terms of gap times^i. 

The following theorem provides sufficient conditions for the existence of a NHIM (the 
proof is given in Appendix |A]) : 

Theorem 1 Consider a < and the region on the n 2 plane bounded by the curve: 



Then this region on n 2 is a two-dimensional (non-isoenergetic) normally hyperbolic invariant 
manifold. 

Note that for a given three dimensional energy surface the NHIM is a (one dimensional) 
trajectory on n 2 . 

C. The existence of a phase space dividing surface having the no-recrossing prop- 
erty 



H 2 = y + 4/3 A* + %A\ . 



(3.12b) 




(3.13) 



where 




(3.14) 



For the Hamiltonian (13.71) we now construct a dividing surface in phase space having 
the "no-recrossing" property. We will describe what this means, as well as the dynamical 
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significance of the dividing surface, in the course of our construction. 

The codimension one non-isoenergetic surface defined by A% = divides the phase 
space into two regions: one associated with the potential well whose minimum is 
(A~i,pi, A~2,p2) = (V — a J 0, 0, 0) and the other associated with the potential well whose 
minimum is (Ai,pi, A 2 ,p2) = (—- \/— a, 0, 0, 0) . The dividing surface restricted to a fixed 
energy surface H = E is given by: 

{—2 —2 ^ 
(A 1 ,p 1 ,A 2 ,p 2 ) \A 1 = 0,H = ^- + ^ + A(3A 2 2 + 8AI = e\ (3.15) 

This dividing surface has two halves: 

{—2 —2 ^ 
(A 1 ,p 1 , A 2 ,p 2 ) \A 1 = 0, H=y + y + 4 ^2 + 8A 4 2 = E, p 1 > j (3.16a) 



and 



{—2 —2 "\ 
(A 1 ,p 1 , A 2 ,p 2 ) \A 1 = 0, H=y + y + 4 ^2 + *Ai = E, p x < j . (3.16b) 

These two halves meet at the NHIM: 

NHIM(E) = | p l9 A 2 ,p 2 ) | A l = 0, # = ^ + 4/3^ + = E, p 1 = o| . (3.17) 

The nature of the NHIM (i.e., the boundary between DS + (E) and DS^(E)) depends on 
both E and f3. The dynamics on the A 2 — p 2 plane is illustrated in Fig. [TJ 

We now argue that DS + (E) and DS_(E) are surfaces having the no (local) re-crossing 
property. These surfaces are defined by A\ = 0. Therefore points on these surfaces leave if 
A 1 ^ 0. We see from flBTTOl that A x = |^ = p x . Therefore on DS + (E) we have Ai > 
and on DS-(E) we have A\ < 0. Trajectories through points on DS + (E) move towards the 
region of phase space associated with the potential well whose minimum is (Ai,pi, A 2 ,p 2 ^ = 
[yf— a, 0, 0, 0) and points on DS_(E) move towards the region of phase space associated 
with the potential well whose minimum is (Ai,pi, ^2,^2) = (—a/— a, 0, 0, 0). 

We denote the directional flux across these hemispheres by <f)+(E) and 4>-(E), respectively, 
and note that 0+(£) + (j).(E) = 0. The magnitude of the flux is \<f>+(E)\ = \<f>-(E)\ = <p(E). 
The magnitude of the flux and related quantities are central to the theory of isomerization 
rates, as discussed in Sec. [TV] 
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IV. PHASE SPACE VOLUMES, GAP TIMES, AND REACTIVE FLUX 



In this section we briefly review the concepts from classical reaction rate theory that will 
be applied to the dynamics of the buckled nanobeam. 

Points in the 4-dimensional system phase space Ai = M 4 are denoted z = 
(pi,p2,Ai, A2) = (p, A) £ M.. The system Hamiltonian is H(z), and the 3 dimensional 
energy shell at energy E, H(z) = E, is denoted S^cM. The corresponding microcanoni- 
cal phase space density is 8(E — H(z)), and the associated density of states for the complete 
energy shell at energy E is 

p(E)= [ dz 5{E - H{z)). (4.1) 

The disjoint regions of phase space separated by DS(i?) are denoted M±; the region 
of phase space corresponding to the potential well whose minimum is ^2,^2) = 

[y/—a, 0, 0, 0) will be denoted by A4+, and that corresponding to the potential well whose 
minimum is (A t ,pi, A 2 ,p2) = {—\/—a 1 0, 0, 0) will be denoted by M... 

The microcanonical density of states for points in region Ai+ is 

P+ (E)= [ dz5(E-H(z)) (4.2) 

JM + 

with a corresponding expression for the density of states P-(E) in .M_. Since the flow is 
everywhere transverse to DS±(£ I ), those phase points in the region M. + that lie on crossing 
trajectories^ 1 ^ (i.e., those trajectories that cross DS±(£')) can be specified uniquely by co- 
ordinates (p, A, ip), where (p, A) £ DS + (-E) is a point on DS + (-E), specified by 2 coordinates 
(p, A), and ip is a time variable. The point z(p,A,ip) is reached by propagating the initial 
condition (p, A) £ DS + (-E) forward for time ifM.^. As all initial conditions on DS + (i?) 
(apart from a set of trajectories of measure zero lying on stable manifolds) will leave the 
region M. + in finite time by crossing DS.(-E'), for each (p, A) £ DS + (i?), we can define the 
gap time s = s(p, A), which is the time it takes for the trajectory to traverse the region A4 + 
before entering the region That is, z(p,A,ip = s(p, A)) £ DS.(-E'). For the phase point 
z(p, A, ip), we therefore have < ip < s(p, A). 

The coordinate transformation z — > (E, ip,p, A) is canonical^ - — , so that the phase space 
volume element is 

d 4 2 = dE dtp da (4.3) 
with da = dpdA an element of 2 dimensional area on the DS(E). 
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The magnitude (f)(E) of the flux through dividing surface DS + (£') at energy E is given 

by 



<f>(E) 



(4.4) 



da 

where the element of area da is precisely the restriction to DS(E) of the appropriate flux 
2-form u corresponding to the Hamiltonian vector field associated with H(z)2£-&t— . The re- 
actant phase space volume occupied by points initiated on the dividing surface with energies 
between E and E + dE is therefore^i^^i^^— 

dE [ da [ dip = dE [ das (4.5a) 
Jds+(e) Jo Jbs+(e) 

= dE(j)(E)s (4.5b) 

where the mean gap time s is defined as 

da s (4.6) 



1 



HE) 



DS+(E) 



and is a function of energy E. The reactant density of states p+(E) associated with crossing 
trajectories only (those trajectories that enter and exit the region A4+— ) is then 

p c + {E) = <f,(E)s (4.7) 

where the superscript C indicates the restriction to crossing trajectories. The result (j4.7p is 
essentially the content of the so-called classical spectral theorem^^^r— . 

If all points in the region Ai + eventually leave that region (that is, all points lie on 
crossing trajectories^ 1 ^) then 

p c + (E) = P+ (E), (4.8) 

so that the crossing density of states is equal to the full reactant phase space density of 
states. Apart from a set of measure zero, all phase points z G M + can be classified 
as either trapped (T) or crossing (C)^. A phase point in the trapped region A^^ never 
crosses the DS(E), so that the associated trajectory does not contribute to the reactive flux. 
Phase points in the crossing region Ai+ do however eventually cross the dividing surface, 
and so lie on trajectories that contribute to the reactive flux. In general, however, as a 
consequence of the existence of trapped trajectories (either trajectories on invariant trapped 
2-torii£*i£ or trajectories asymptotic to other invariant objects of zero measure), we have the 
inequality^^^ 

p°(E) < p + (E). (4.9) 
14 



If p+(E) < p + (E), then it is in principle necessary to introduce corrections to statis- 



21 



48 



crossing and reactive 
and results for the 
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Note that, if the strict inequality 



tical estimates of reaction rates^i&^Sr— . Numerical computa tion o 
densities of states for the HCN molecule are discussed in refs 
Hamiltonian isokinetic thermostat are discussed in ref. 
p%(E) < p+(E) holds, then the system dynamics cannot be ergodic on the energy shell at 
energy E. The equality p+(E) = p+(E) is therefore a necessary condition for ergodicity, one 
that can be checked numerically. 



A. Gap time and reactant lifetime distributions 

The gap time distribution, V(s; E) is of central interest in unimolecular kinetics^ 1 ^: the 
probability that a phase point on DS+(i£) at energy E has a gap time between s and s + ds 
is equal to V(s; E)ds. An important idealized gap distribution is the random, exponential 
distribution 

V(s;E) = k(E)e~ k{E)s (4.10) 

characterized by a single decay constant k (where k depends on energy E), with correspond- 
ing mean gap time s = k^ 1 . An exponential distribution of gap times is usually taken to be 
a necessary condition for 'statistical' behavior in unimolecular reactions^ 1 ^— . 

The lifetime (time to cross the dividing surface DS.(-E')) of phase point z(p,A,ifj) is 
t = s(p,A) — ip, and the corresponding (normalized) reactant lifetime distribution function 
P(t; E) at energy E i^Stm^. 

P(t; E) = -— Prob(t > t'; E) ^ ^ (4.11a) 

+oo 

dsV(s;E) (4.11b) 



dt 

1 



s 

where the fraction of interesting (reactive) phase points having lifetimes between t and t + dt 
is P(i; E)dt. It is often useful to work with the unnormalized lifetime distribution F, where 
F{t- E) = sP(t; E). 

Equation ( 14.1 lap gives the general relation between the lifetime distribution and the frac- 
tion of trajectories having lifetimes greater than a certain value for arbitrary ensembles^— . 
Note that an exponential gap distribution (14 . 1 f) implies that the reactant lifetime distribu- 
tion P(i; E) is also exponent ia^^ 1 ^^— ; both gap and lifetime distributions for realistic 
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molecular potentials have been of great interest since the earliest days of trajectory simu- 
lations of unimolecular decay, and many examples of non-exponential lifetime distributions 
have been found^^r— . 



B. Reaction rates and the inverse gap time 

The quantity 



k RRKM (E) = (4.12) 



is the statistical (RRKM) microcanonical rate for the forward reaction (trajectories crossing 
DS + ) at energy E, the ratio of the magnitude of the flux 4>(E) through DS + (i?) to the total 
reactant density of states^^. 
Clearly, if p + (E) = p9(E), then 



k RRKM (E) = = (4.13) 



1 



7 



s 



the inverse mean gap time. In general, the inverse of the mean gap time is 

*4 = f) (4,4a) 



lRRKM 
K f 



P + (E) 

Lp5(^)J 



(4.14b) 



> k RRKM . (4.14c) 

The inverse gap time can then be interpreted as the statistical unimolecular reaction rate 
corrected for the volume of trapped trajectories in the reactant phase spaced 



C. Reactive flux correlation function 



The discussion of reactive fluxes across the phase space dividing surface separating re- 
actant from product and of gap times provides a theoretical framework for analyzing the 
lifetime distribution of an ensemble of trajectories initiated in the reactant well at constant 
energy, where the lifetime refers to the time to the first crossing of the dividing surface. 

Another approach to isomerization kinetics considers an equilibrium (canonical or micro- 
canonical) ensemble of reactants and products. The regression hypothesis relates the total 
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relaxation rate for an initial perturbation of the equilibrium populations to the autocorre- 
lation function of spontaneous population fluctuations^. Standard analysis^ then provides 
a relation between the isomerization rate, when the latter exists, and the computationally 
tractable quantity /C(t) given in terms of the reactive flux across the barrier: 

m = -i-(g(0)<%(0) - g*]e + [g(f)]> ~ < -^. (4.15) 

In this expression, q = Ax, the reaction coordinate, and the DS is determined by symmetry, 
so that the critical value = 0. We also have 

- = k f + k b = 2k f . (4.16) 

In our calculations the ensemble average (■ ■ •) corresponds to a microcanonical ensemble 
(average over the entire energy shell 0+ is the characteristic function for the configura- 
tion space region A 2 > 0, and equilibrium fractions £1X6 — X = 1/2. In the limit t -> 0, 

the right hand side of equation (I4.15P is just twice the statistical rate /^ RKM , eq. (I4.12p . 

Operationally, in principle we must sample the DS without regard to the sign of the initial 
velocity q(0). A trajectory contributes to the average ( 14. 15}) at time t: 

(i) Only if the phase point is in the product well (q > 0) at time t, 

(ii) With a sign (±) determined by the initial sign of q. 

The right hand side of (14. 15j) decays to zero as t — > oo, as trajectories initially crossing 
from product to reactant (q < 0) eventually return to the product side, leading to can- 
cellation. If the right hand side of ( 14.151) exhibits a so-called 'plateau' region in which it 
is approximately constant, followed by exponential decay, then an isomerization rate can 
be extracted from the computation. This behavior indicates a well-defined separation of 
timescales: trajectories remain trapped in either well for long times with only infrequent 
transitions (crossing of the DS) between wells. On the other hand, if the reactive flux 
correlation functions exhibits oscillatory decay, then no rate constant exists at the energy 
in question. Pioneering computations of flux correlation functions for a number of 2 DoF 
dynamical models for isomerization were made by DeLeon and Berne^ 1 ^. 

In practice, we exploit the symmetry of the potential, and sample only initial conditions 
with Ai > 0. If the fraction of the phase points in the product well A\ > at time t is W(t), 
then the fraction of the phase points in the reactant well, A\ < 0, at time t is 1 — W(t). 
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As the potential is symmetric about A\ = 0, reversing the initial sign of q leads to a set of 
symmetry-related trajectories with the occupancy of the two wells now 1 — W(t) and W(t), 
respectively Adding the contributions of trajectories to (14.151) with appropriate sign yields 
a result proportional to 2W(t) — 1, which can be calculated from W(t) directly. 

A connection between the gap time and the reactive flux approaches to isomerization 
kinetics was established by Straub and Berne in their work on the "absorbing boundary" 
method for computing isomerization rates^^. Assuming that there are no correlations 
between successive crossings of the DS for a given trajectory (i.e., assuming "chaotic" dy- 
namics), then the single-passage gap/lifetime distribution can be used to derive an expression 
for the reactive flux, and hence the associated isomerization rate^. 
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V. PARAMETER VALUES AND ENERGY SCALES 



A. Physical parameters 

We study a 2-mode truncation of the dynamics of a silicon nanobeam having rectangular 
cross section under compressive stress, subject to hinged boundary conditions. The following 
physical parameter values are used 9 - - — : elastic modulus Q = 1.3 x 10 11 J/m 3 ; uncompressed 
length L = 5 x 10~ 8 m; width w = 2 x 1CT 9 m; depth d — 1 x 1CT 9 m; density p = 2330 
kg/m 3 . For these parameters, the critical value of the strain e c , obtained by solving the 
implicit equation fl2T4bl . is e c = -0.000329. 

B. Strain values: 3 cases 

For a beam described by the physical parameters listed above, we consider the dynamics 
for 3 values of the compressive stress. The corresponding strain values and associated param- 
eters are given in Table HT1 Strain values for the three cases are: case I, e = —0.00065840 ~ 
2 x e c ; case II, e = -0.00197520 ~ 6 x e c ; case III, e = -0.001419692 ~ 4 x e c . The strain 
value for case III is chosen so that the energy of the index 2 saddle lies just above the pair 
of index 1 saddles at = —0.0001. 

Contour plots of the potential (I2.19p for the 3 cases considered are given in Figure EJ 
Setting coordinate A 2 = in potential f)2.19p . we obtain a bistable potential which is a 
function of the 'reaction coordinate' A\. In Table HT1 we give the value of the barrier height 
AE for each of the 3 cases (degrees K). It can be seen that the barrier heights are compa- 
rable to thermal energies ~ 100 K for all cases. We have also estimated the magnitude of 
vibrational quanta hu associated with oscillations of the beam along the reaction coordinate 
at the potential energy minimum; these energies are given in Table [III We have hco / k^T 1 
for T > 100 K. 
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VI. RESULTS AND DISCUSSION 



A. Reactive flux, phase space volumes and ergodicity 

We have computed reactive fluxes associated with the symmetry-determined DS A 2 = 
for each of the 3 cases listed in Table HT1 at 3 energies: E = 10~ 9 , E = 1CT 8 and E = 1CT 7 . 
For cases II and III, where the coordinate origin (0, 0) is an index 2 saddle on the potential 
energy surface flanked by a pair of index 1 saddles, we have also performed computations 
at several values of E < 0. For case II, where the energy of the index 2 saddle is well 
above the pair of index 1 saddles, we have used 3 additional energies: E = —2.12 x 10 -7 , 
E = —2 x 10~ 7 and E — — 1 x 10 -7 . For case III where the energy of the index 2 saddle is 
only just above the pair of index 1 saddles we consider 2 additional energies: E = —4 x 10 -9 
and E = —2.5 x 10~ 9 . In each case the lowest energy is close to the index 1 saddle energy. 
Our numerical results, obtained via Monte Carlo sampling of the DS, are presented in Tables 
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IHIfTVl For further details on the numerical methods used in these computations, see ref 

Numerical results for mean gap time s, reactive flux (p + (E), reactant volume p%(E) = 
sx0 + , reactant density of states p + (E), pulse decay constant k (see below) and the statistical 
isomerization rate /c^- RKM for each case and energy studied are given in Tables IIIlMVl 

For the simple quartic potential obtained by setting Ai = 0, action integrals (fluxes) 
12(E) for motion in the invariant plane Il 2 can be computed explicitly as a function of total 
energy E in terms of complete Elliptic integrals. These analytical expressions, not reported 
here, have been used as a check on our numerical calculations. 

Our results show that p+(E) < p+(E) in all instances; in the majority of cases the phase 
space volume swept out by reacting (crossing) trajectories is considerably smaller than the 
full classical density of states associated with the reactant region of phase space. This means 
that, for the stress values and energies studied here, the buckled nanobeam dynamics is very 
far from being ergodic. Ergodicity is usually taken to be a necessary (but by no means 
sufficient) condition for the applicability of statistical theories of reaction rates. 

Some representative trajectories for the 3 cases are shown in Figure [3j For each 
case/energy we present two plots: one shows 20 trajectories initiated on the DS and fol- 
lowed until the first recrossing of the DS, while the other shows a single trajectory followed 
for 200 crossings of the DS. It is clear by inspection of the single trajectory plots that the 
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dynamics is far from ergodic on the timescale considered; the trajectories appear to be 
quasiperiodic or weakly chaotic at most. 

B. Gap time distribution 

Both gap time distributions V(t) and associated (unnormalized) lifetime distributions 
F(t) have been computed for all cases. The functions V(t) and log[F(t)] are plotted in 
Figure H] for the lowest energy in each case. The results shown represent the range of 
behavior found for the various cases and energies we have studied. 

For case I, E — 10~ 7 , the gap time distribution essentially consists of a single 'pulse' 
associated with trajectories that enter the well, exhibit a single turning point in the reaction 
coordinate, and then exit through the dividing surface DS_. The smallest gap time is 
nonzero, reflecting a delay corresponding to the time it takes for a point on the shortest lived 
trajectory to reach the turning point and then return to the DS. The lifetime distribution of 
the set of trajectories in the pulse is however well described by a single exponential decay, 
with a decay rate (denoted k) that is much faster than either the RRKM rate ^ RKM or the 
inverse of the mean gap time. 

For cases II and III, the structure of the gap time distribution is more complex, consisting 
of multiple pulses. Several early pulses have comparable amplitudes, with amplitude not 
necessarily decreasing monotonically with gap time. In such cases the computation of decay 
rates for individual pulses is less accurate due to the possibility of overlapping pulses. At 
long gap times the distribution becames smooth, with typically nonexponential decay. 

Representative pulse decay constants k are listed in Tables IllhWl In all instances k ^> 

lRRKM 
K f 

We note that the gap time distributions seen here are reminiscent of the 'epistrophic' 
patterns of ionization times seen in the work of Mitchell and Delos 7 ^. 

C. Reactive flux correlation function 

The quantity fC{t) (eq. (I4.15P ) is plotted in Figure for three different cases/energies. 
The behavior seen in these plots is again typical of all the cases we have examined: JC(t) 
exhibits oscillatory decay over a much longer timescale than the gap time decay constant k. 
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The absence of a 'plateau' region means that an isomerization rate constant does not exist 
for our nanobeam model in the physical regime studied. 

D. Gap time distribution on the DS 

To further explore the isomerization dynamics of the nanobeam, we examine the distri- 
bution of gap times on the dividing surface. That is, we plot contours of the gap time s as 
a function of coordinates (A 2 ,p 2 ) on DS + . 

A set of representative plot is shown in Figure El corresponding to case I, E — 10 -7 , 
case II, E = 10 -9 and case III, E = 10~ 9 . Note that the contour plots are invariant with 
respect to the inversion operation (A 2 ,p 2 ) (~A 2 , — p 2 ). For each case the gap time s is 
also plotted along the line p 2 = 0. 

A significant finding is that for case I the gap time is a smooth function of location on the 
DS; although there appear to be no singularities inside the boundary (NHIM) with divergent 
gap times, closer examination (not shown here) confirms the existence of a small region near 
the NHIM (the boundary of the DS) where gap times are longer, presumably associated 
with second (and later) pulses that do not cross the TS on first approach but bounce back 
again one (or more) times. The absence of 'fractal' patterns such as those seen in previous 
studies^ indicates that the intra-well dynamics is extremely simple for case I. 

The gap time contours on the DS for cases II and III show a more typical^ fractal 
arrangement. Initiial condition on the DS for which the gap time diverges presumably lie 
on the stable manifold of either the NHIM or of a periodic orbit confined to the reactant 
well. By symmetry, trajectories on the line p 2 = having divergent gap times must lie on 
both the stable manifold (in forward time) and unstable manifold (in backward time) of the 
NHIM, and hence lie on homoclinic orbits. 
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VII. CONCLUSIONS AND OUTLOOK 



In this paper we have studied the classical (Euler-Bernoulli) mechanics of a 2-mode 
truncation of the dynamics for a buckled nanobeam with rectangular cross section subject 
to compressive stress. The physical parameters used correspond to a Silicon nanobeam^r— . 
In the stress regime studied, the first transverse displacement mode has become unstable, 
while the second mode is either stable or unstable, depending on the value of the strain. The 
resulting beam Hamiltonian has the same form as model 2 DoF systems previously studied 
in chemistry, which describe a bistable reaction (isomerization) coordinate coupled to an 
additional transverse or 'bath' mode. 

We have applied methods from reaction rate theory to characterize the nanobeam 'iso- 
merization' dynamics. For the beam model considered, the dividing surface separating 
'reactant' and 'product' configurations for the buckled beam is completely determined by 
symmetry (coordinate A\ — 0). Using exponential dichotomies, we have proved that, for a 
specified range of the energy, the boundary of the associated dividing surface is a Normally 
Hyperbolic Invariant Manifold (NHIM). 

We have computed reactive fluxes, mean gap times and reactant phase space volumes 
for 5 stress values at several different energies. In all cases the phase space volume swept 
out by crossing trajectories is considerably less that the reactant density of states, proving 
that the dynamics is highly nonergodic. The associated gap time distributions consist of 
single 'pulses' of trajectories. Computation of the reactive flux correlation function shows 
no sign of a plateau region; rather, the flux exhibits oscillatory decay, indicating that, for 
the 2-mode model in the physical regime considered, a rate constant for isomerization does 
not exist. 

Problems for future work include study of the dynamical influence of additional 'bath' 
modes, and the investigation of quantum effects at low temperatures. 
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Appendix A: Proof of the Existence of a Normally Hyperbolic Invariant Manifold 



Roughly speaking, a normally hyperbolic invariant manifold has the property that, under 
the dynamics linearized about the invariant manifold, growth rates in directions transverse 
to the invariant manifold dominate the growth rates of directions tangent to the invariant 
manifold. (For some background on normally hyperbolic invariant manifolds see refs 
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An account of Fenichel's approach to the theory, as well as some history and examples can 
be found in ref. |28|. ) The dynamics on II2 are completely integrable. For (3 > it consists 
entirely of periodic orbits, and the growth rates (e.g. Lyapunov exponents) associated with 
all orbits are all zero. For (3 < it consists entirely of periodic orbits, except for the saddle 
point at the origin connected by a pair of homoclinic orbits. The periodic orbits all have zero 
growth rates, and the saddle point and the homoclinic orbits will be discussed separately. 
We will show that the growth rates transverse to II 2 are exponential. Hence they dominate 
the growth rates tangent to II2. 

We linearize ( 13 . 1 0[) about IT 2 and evaluate the resulting equations on an arbitrary tra- 
jectory on n 2 . Since n 2 is a plane, and is described in a global coordinate system, lin- 
earization about Il 2 is particularly easy. Coordinates describing the directions normal to n 2 
are (Ai,pi), and Il 2 is defined by Ai = pi = in the four dimensional phase space with 
coordinates (A 1 ,pi, A 2 ,p2)- 

Therefore, the linearized dynamics about n 2 is obtained by retaining terms only linear 
in (Ai,pi): 

Ai=px, (Ala) 
pi = -2Ai (a + AAl) , (Alb) 
A 2 = p 2 , (Ale) 
P2 = -8A 2 {(3 + AAj) . (Aid) 

The linearized dynamics normal to Il 2 are given by: 

Ai=Pi, (A2a) 
pi = -2Ai (a + AA\{t)) . (A2b) 

where A 2 (t) is a component of a trajectory on Il 2 , i.e. it is the A 2 component of a trajectory 
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of: 



A 2 =p 2 , (A3a) 
p 2 = -8A 2 (p + AAj) . (A3b) 



We now will show that the linear, nonautonomous equation (IA2[) has two linearly inde- 
pendent solutions; one exponentially growing in time, and the other exponentially decaying 
in time. This is accomplished by showing that (1A2j) has an exponential dichotomy, and these 
conditions will depend on a, (3, A^it). 

The numbers that are typically used to quantify (linearized) growth rates of trajectories 
are the Lyapunov exponents. However, for nonautonomous ordinary differential equations 
exponential dichotomies are often more convenient and facilitate the proof of certain results. 
A discussion of the relationship between Lyapunov exponents and exponential dichotomies 



can be found in, e.g., refs 
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a. Exponential Dichotomy First, we will provide some background on the notion of ex- 
ponential dichotomy that is particular to our situation, i.e. two dimensional time-dependent 
ordinary differential equations. The standard reference on exponential dichotomies is ref. 
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Consider the linear ordinary differential equation with time dependent coefficients 

x = L(t)x, (A4) 

where x = (x, y) and the 2x2 matrix L(t) is a continuous function of t, and suppose 
X(t) is the fundamental solution matrix of (1A4I) . Let || ■ || denote a matrix norm, such as 
the maximum of the absolute values of the matrix elements. Then (IA4I) is said to possess 
an exponential dichotomy if there exists a rank-one projection operator P and constants 
Kx, K 2 , ai, a 2 > 0, such that 

|| X{t)PX-\r) || < K 1 exp(-o;i(t-T)), t > r, (A5a) 
|| X(t) (id- P)X~ 1 (t) (I < K 2 exp(a 2 (t-r)) , t < r. (A5b) 

The condition that the projection operator P has rank one means that, of the two lin- 
early independent solutions of (1A4I) . one is exponentially growing and one is exponentially 
decaying. 
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Verifying that (1A2|) has an exponential dichotomy requires us to solve for the fundamental 
solution matrix X(t). In general, this is not possible. Instead, we will use results on 
"roughness of exponential dichotomies"— 1 ^. The relevant result is as follows. Suppose flA4j) 
has the form: 

±={A(t)+B{t))x, (A6) 

and suppose that the equation: 

x = A(t)x, (A7) 
has an exponential dichotomy with constants K\, K 2 , oci, a 2 - If 



, . , .. K x K 2 . 
sup || B(t) || — + — < 1, 



(A8) 



then ( 1A6I) has an exponential dichotomy with constants K[ = K 2 = K3 > and a[ = a 2 
a 3 > 0. 

We now apply these ideas to (1A2j) . Rewriting this equation in the form of (1A4j) gives: 



Pi 



1 

-2a 





L4|(t) 



Pi 



(A9) 



We introduce a linear change of coordinates that diagonalizes the first matrix in this expres- 
sion: 



Mil 






M 




U 1 







where 



1 



1 



T -i 



Then (IA9D becomes: 







( 









-v^2^ -1 



2 v /r 2o : I _ v /32^ 1 



441(0 



(A10) 



(All) 



V^2a 1 

-v^^ / V^2a~ I -4^2(t) -4A|(t) 



.r 
2/ 



First we consider: 



1 

The fundamental solution matrix is given by: 

A/ = 2at q 



)l 






W 



(A12) 



(A13) 



X{t) 



e" 
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-2a t 



(A14) 




We take as the projection matrix: 

p = i : ; i (Ms) 

Then, using (lAISp and (1A14[) . we have: 

X{t)PX-\r) = ° ° 

I e -V=2H(t-r) 

and 

AT (t) (id -P)X^(r) = I (A16b) 




It follows that ( 1A13I) has an exponential dichotomy with K\ = K2 = 1 and a\ = 02 



y/—2a. We now need to check (1A8|) . For (1A12[) this condition takes the form: 

4 sup teR |A|(t)| < j ^ Al7 ^ 
—a 

The quantity sup tgR is the maximum value that A\(t) attains along a trajectory, 

and this can be precisely computed, as we now show. 

In Section IIIIBI we showed the the dynamics in the A 2 — P2 plane is Hamiltonian, with 
Hamiltonian given by: 

H 2 = ^ + + 8Ai (A18) 
Trajectories lie on level curves of the Hamiltonian B.2- 

A + A(5A\ + 8AI = E, (A19) 

where 

E > for > 0, 



(A20) 

E > -\I3 2 for < 0. 
The quantity sup f6K |A|(t)| corresponds to the largest "turning point" of a trajectory, 
i.e., the largest value of A\ that intersects the A 2 axis. An equation for this quantity can be 
obtained by setting p2 = in ( 1A19I) . Doing so, and rearranging terms, gives the following 
quadratic equation for A%: 

At + -A\-- = 0. (A21) 
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The solution of this equation is: 



a I = -j±\Vp*+ze ( A22 ) 



and we will take the plus sign in front of the square root since we are seeking the largest 
root. It is also useful to note that this expression is an increasing function of E since: 

dA 2 1 /A . 

> 0. (A23) 



dE 4 v//3 2 + 2E 

Therefore we seek the largest value of E such that (IA17I) is satisfied. This is obtained by 
equating the expression for A 2 , = — | + \\/ ft 2 + 2E to — j, which gives: 

-P+ v 7 /? 2 + 2E = -a. (A24) 

and then solving this expression for E: 

E max = ^ ( 1 - 2^V (A25) 
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Equilibrium point 


Eigenvalues 


Energy of the equilibrium 


{A 1 ,p 1 ,A 2 ,p 2 ) = (0,0,0,0) 


±V-2a, ±2^-2(1 





(A 1 ,p 1 ,A 2 ,p2) = {±V^a, 0,0,0) 


±2^0, ±2y/2(a-P) 




{A 1 ,p 1 ,A 2 ,P2) = (o,o,±^,o) 


±Vft ±V 2 (0-«) 


-¥ 2 



TABLE I: The location of the equilibria, the eigenvalues of the matrix associated with the lin- 
earization of Hamilton's equations about the equilibria, and the (total) energy of the equilibrium. 
In the last two rows both equilbrium points have the same four eigenvalues. 



Case 


e 


e 


a 


/3 


AE/k B [K] 


hoo/k B [K] 


I 


-0.00065840 


-0.00032942 


-0.00032898 


0.00065928 


50.98 


0.092 


II 


-0.00197520 


-0.00033029 


-0.00164491 


-0.00065404 


1274.4 


0.206 


III 


-0.00141969 


-0.00032992 


-0.00108977 


-0.00010000 


559.4 


0.168 



TABLE II: Strain values e and associated parameter values e, a and /3 for the 3 cases used in our 
computations. Also shown are barrier heights AE and estimates of the size of vibrational quanta 
for beam oscillations about the energy minimum (degrees K). 



Energy 


s 


ME) 




P+(E) 


k = s- 1 


fMRKM 


K 


le-09 


357.460 


0.00000009 


0.000031 


0.00085 


0.00280 


0.0001013 


0.0254 


le-08 


266.663 


0.00000086 


0.000230 


0.00103 


0.00375 


0.0008351 


0.0248 


le-07 


174.982 


0.00000832 


0.001456 


0.00210 


0.00571 


0.0039613 


0.0224 



TABLE III: Computational results for case I: a = -0.00032898, /3 = 0.00065928. For discussion 
see Sec. EH 
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Energy 


s 


ME) 


P° + {E) 


P+(E) 


k = s~ 1 


jLRRKM 

f 


K 


-2.12e-07 


7729.396 


0.00000023 


0.001789 


0.00897 


0.00013 


0.0000258 


0.0645 


-2e-07 


1258.603 


0.00000172 


0.002160 


0.00914 


0.00079 


0.0001878 


0.0406 


-le-07 


385.348 


0.00001488 


0.005734 


0.01034 


0.00260 


0.0014386 


0.0481 


le-09 


166.588 


0.00003186 


0.005308 


0.01137 


0.00600 


0.0028028 


0.0312 


le-08 


165.124 


0.00003384 


0.005588 


0.01132 


0.00606 


0.0029886 


0.0454 


le-07 


141.782 


0.00004771 


0.006764 


0.01183 


0.00705 


0.0040325 


0.0492 



TABLE IV: Computational results for case II: a = -0.00164491, f3 = -0.00065404. For discussion 
see Sec. IVII 



Energy 


s 


ME) 


P C + (E) 


P+(E) 


k = s- 1 


ilRRKM 


K 


-4e-09 


4076.711 


0.00000032 


0.001308 


0.00574 


0.00025 


0.0000559 


0.0618 


-2.5e-09 


1567.930 


0.00000083 


0.001304 


0.00579 


0.00064 


0.0001436 


0.0427 


le-09 


502.955 


0.00000236 


0.001187 


0.00590 


0.00199 


0.0004001 


0.0402 


le-08 


307.775 


0.00000497 


0.001528 


0.00600 


0.00325 


0.0008280 


0.0621 


le-07 


195.281 


0.00001941 


0.003791 


0.00680 


0.00512 


0.0028533 


0.0410 



TABLE V: Computational results for case III: a = -0.00108977, f3 = -0.0001. For discussion see 
Sec. El 
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Figure captions 



FIG. 1: Phase space portraits in the invariant A 2 — p 2 plane, (a) /3 > (b) < 0. 
FIG. 2: Contour plots of the 2- mode nanobeam potential eq. (I3.8P for the 3 compressive 
stress values considered in this paper. The contour values shown include the particular 
energies at which the dynamics was studied, (a) a = —0.00032898, (3 = 0.00065928; (b) 
a = -0.00164491, (3 = -0.00065404; (c) a = -0.00108977, (3 = -0.0001. 

FIG. 3: Plots of trajectories initiated on the DS. (a), (b) Case I, energy E = 10 -7 ; (c), (d) 
case II, energy E = 10~ 9 ; (e), (f) case III, energy E = 10~ 9 . Panels (a), (c) and (e) each 
show 20 trajectories followed until the first recrossing of the DS, while panels (b), (d) and 
(f) show single trajectories followed for 200 crossings of the DS. 

FIG. 4: Gap time distribution V(t) and the logarithm of the associated (unnormalized) 
lifetime distribution E(t). (a), (b) case I, energy E = 10~ 7 ; (c), (d) case II, energy E = 10~ 9 ; 
(e), (f) case III, energy E = 10 -9 . 

FIG. 5: Reactive flux correlation function JC(t) versus t. (a) case I, energy E = 10 -7 ; (b) 
case II, energy E = 10 -9 ; (c) case III, energy E = 10~ 9 . 

FIG. 6: Contours of the gap time s as a function of coordinates (^2,^2) on the dividing 
surface DS + (panels (a), (c), (e)) and as a function of A 2 along the line p 2 = (panels (b), 
(d) and (e)). (a), (b) case I, energy E = 10~ 7 ; (c), (d) case II, energy E = 10~ 9 ; (e), (f) 
case III, energy E = 10~ 9 . 
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FIGURE 2 
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FIGURE 4 
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